*****************************************************************
* Replication directory for                                   ***
* Prime locations                                             ***
* by Gabriel M. Ahlfeldt, Thilo N.H. Albers, Kristian Behrens ***
* Published in American Economic Review: Insights             ***
*****************************************************************
* 01/2025
* Stata
version 17.0


* This do file generates maps of prime locations by Global City

* Create directory
	capture mkdir "$figures_App/GlobalCities"
	capture mkdir "$figures_App/GlobalCities/GISMAPS"

* Load grid output
	u "$temp/125CitiesPLs/grid125_PL_output.dta", clear	

* Merge ranks
	merge m:1 PLID metro_id using "$temp/PLID125_metro_rank"
	drop _m 

* Merge metro name 
	merge m:1 metro_id using "$data_125cities/METRO_LEVEL_COVARIATES/metro_data.dta", keepusing(metro)
	ren metro metro_name
	drop if _m == 2
	drop _m
* Merge population
	merge 1:1 grid_x grid_y metro_id using "$data_125cities/POP/GRID_POP.dta", keepusing(popdens)
	drop if _m == 2
	drop _m
	replace popdens = 0 if popdens == .
* Merge photos	
	merge m:1 grid_x grid_y  using "$data_125cities/PHOTOS/PHOTOCOUNT_250M.dta"
	drop if _m == 2
	drop _m
	ren count photos
	replace photos = 0 if photos == .
	replace photos = 0 if developable == 0

* Create PL and save for mapping later
	preserve 
	keep if PL == 1
	save "$temp/PL125_grids", replace
	restore 	
	
	preserve 
	keep if developable == 0
	keep cell_id developable metro_id
	save "$temp/undev125_grids", replace
	restore 		
	
* Define metro list
	sum metro_id 
	global MetroMin = r(min)
	global MetroMax = r(max)
	
* Save data	
	capture mkdir "$temp/GIS125_using"
	save "$temp/GIS125_using/base", replace

* Load metro list	
	qui import delimited "$data_125cities/METRO_LEVEL_COVARIATES/metrolist125.csv", clear
	qui tab metro_id
	local Nmetro = r(N) // Total number of cities in list
	local count = 1
	levelsof  metro_id, local(METROIDS)		
	 
* Loop through metros
	foreach num of local METROIDS { // loop over metro	
		display "...generate spatial data for metro `num'..."
		qui shp2dta using "$data_125cities/GIS/CELLS/CELL_`num'.shp", database("$temp/GIS125_using/db_using_`num'") coordinates("$temp/GIS125_using/coord_using_`num'") replace // Full grid with employment data
	}	
* Use the coordinates shape and restrict it to the PL grid cells and save for each city a file	
	foreach num of local METROIDS { // loop over metro	
			preserve
			display "...Preparing PLs of metro `num'..."
			qui u "$temp/GIS125_using/coord_using_`num'", clear 
			qui merge m:1 _ID using  "$temp/GIS125_using/db_using_`num'", keepusing(cell_id) // merge ID variable for mapping
				qui drop _m
			qui merge m:1 cell_id using  "$temp/PL125_grids" // merge ID variable for mapping	
				qui keep if _m == 3	
				qui drop _m
				qui keep if metro_id == `num'
				qui save "$temp/GIS125_using/PL_coord`num'.dta", replace
			restore
		}	
* Plot maps		
foreach num of local METROIDS { // loop over metro	
	* Preserve data
	* Merge ID from generated database	
		qui u "$temp/GIS125_using/base", clear
		qui keep if metro_id == `num'
	* Check for ensuring correct scale
		qui sum grid_x
			local xrange = r(max)-r(min)
		* display `xrange'
		local width = int(`xrange'  /40000*2000)
		qui merge 1:1 cell_id using  "$temp/GIS125_using/db_using_`num'", keepusing(_ID) // merge ID variable for mapping
			qui drop if _m == 2
			qui drop _m
	* Now plot maps 
		display "...plotting metro map `num'..."
		global PCT -1.5 -0.5 0.5 1.5 5 50 100 1000 1000000 // defines the cutoffs in the legend
			qui gen labX = PLX+1500
			qui gen labY = PLY // 1500
			qui replace photos = -1 if developable == 0 
			* qui keep if developable == 1
			local title = metro_name[1]
			spmap photos using "$temp/GIS125_using/coord_using_`num'", id(_ID) fcolor(gs14 gs12 gs10 gs8 gs6 gs5 gs4 gs2 ) ocolor(white ..) osize(0.00000 ..)   clmethod(custom)  clbreaks($PCT)  	 legend(off)  /// title("`title'", size(small)) legtitle("Employment")		legend(off) 	/// Purples  clmethod(custom)  clbreaks($PCT)  clmethod(q) cln(5) 
			/// point(xcoord(PLX) ycoord(PLY) size(vsmall) shape(x) ) ///
			label(xcoord(labX) ycoord(labY) label(PLrankByCity)   color(red) size(4) ) ///
			legend(size(3) region(fcolor(gs16))) ///
			polygon(data("$temp/GIS125_using/PL_coord`num'") fcolor(maroon%75) ocolor(red)  osize(0.00000 ..) ) 
* Save maps for toolkit
			qui graph export 	"$figures_App/GlobalCities/GISMAPS/GISMAP_PL_metro`num'.png",  height(1000)  width(1000)  replace // height(1200) width(`width') 
	}
	
* Script ends